CFA como modelo mixto

Ajustando un CFA para obtener resultados de un MLM

Author

https://dacarras.github.io/

Modified

November 14, 2024

Resumen

  • Simulamos datos multinivel con componentes de varianza de \(\sigma^2_{w} = .1\) y \(\sigma^2_{b} = .9\), y con media cero. Generamos datos para 500 clusters, con 5 observaciones para cada cluster, dando un total de 25000 observaciones. Estos datos, se pueden bajar del siguiente link.

  • Cada cluster lo simbolizamos con la letra j, y a cada observación con la letra i.

  • Los datos generados se encuentran en formato stacked o long, donde cada fila es una observación, indexada en un cluster. Estos datos los volteamos, para producir una representación de las mismas observaciones en formato wide o hacia el lado. En este caso, se emplean 5 columnas (i1-i5) para alojar los valores de cada observación por cada cluter j.

  • Empleando a la libreria lavaan, ilustramos como ajustar un modelo MLM sobre los datos en stacked, y ajustamos un modelo CFA sobre los datos en wide, recuperando los mismos resultados en cada caso, hasta a dos decimales.

Flujo de trabajo

    1. Datos simulados
    1. Preparar los datos
    1. Ajustar Modelos
    1. Tabla comparada de resultados

1. Datos Simulados

\[y_{ij} = \gamma_{00} + \mu_{0j} + \epsilon_{ij}\]

Donde,

  • \(y_{ij}\) = respuesta \(y\) para la observación \(_{i}\), en el cluster \(_{j}\)
  • \(\gamma_{00}\) = media entre todos los clusters
  • \(\mu_{0j}\) = desviaciones de las medias de los clusters, a la gran media \(\gamma_{00}\)
  • \(\epsilon_{ij}\) = desviación de las observaciones con respecto a la media de su propio cluster
  • \(\mu_{0j} \sim \mathcal{N}(0,\,\sigma^{2}_{b})\) = donde \(\mu_{0j}\) es variable de media cero, y varianza \(\sigma^{2}_{b}\)
  • \(\sigma^{2}_{b}\) = .9
  • \(\epsilon_{0j} \sim \mathcal{N}(0,\,\sigma^{2}_{w})\) = donde \(\epsilon_{ij}\) es variable de media cero, y varianza \(\sigma^{2}_{w}\)
  • \(\sigma^{2}_{w}\) = .1
  • \(_{j}\) = 500 clusters
  • \(_{i}\) = 5 observaciones para cada cluster

Datos simulados:

Nota

Los datos simulados fueron generados empleando al software Mplus (ver Anexo).

2. Preparar datos

Abrimos los datos simulados, y generamos los datos en formato wide

# ubicación de los datos en url
url_file_long <- url('https://github.com/dacarras/edu4046_cfa/raw/main/data_long.rds')

# datos en formato stacked
data_long <- readRDS(url_file_long)

# datos en formato wide
data_wide <- data_long %>%
             mutate(id_ij = paste0('i',rep(seq(1:5),500))) %>%
             tidyr::pivot_wider(
              id_cols = -id_i,
              names_from = 'id_ij', 
              values_from = 'y') %>%
             dplyr::glimpse()
Rows: 500
Columns: 6
$ id_j <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19…
$ i1   <dbl> -0.649609, 1.810764, -1.365712, -0.555940, -1.024259, -0.441473, …
$ i2   <dbl> -0.513657, 0.345480, -1.018105, -0.355253, -1.083364, -0.886687, …
$ i3   <dbl> -0.697763, 0.818963, -1.361241, -0.655227, -1.675723, -0.422175, …
$ i4   <dbl> 0.009125, 0.982486, -1.054176, 0.053315, -1.880564, 0.001164, 2.1…
$ i5   <dbl> 0.344252, 0.929602, -0.563126, -0.030324, -1.281889, -1.597010, 1…

Cómo se ven los datos generados, y los datos en formato wide, para los 4 primeros clusters.

# datos en formato stacked
data_long %>%
dplyr::slice(1:20) %>%
knitr::kable(., digits = 4)
id_i id_j id_ij y
1 1 1 -0.6496
2 1 2 -0.5137
3 1 3 -0.6978
4 1 4 0.0091
5 1 5 0.3443
6 2 1 1.8108
7 2 2 0.3455
8 2 3 0.8190
9 2 4 0.9825
10 2 5 0.9296
11 3 1 -1.3657
12 3 2 -1.0181
13 3 3 -1.3612
14 3 4 -1.0542
15 3 5 -0.5631
16 4 1 -0.5559
17 4 2 -0.3553
18 4 3 -0.6552
19 4 4 0.0533
20 4 5 -0.0303
# datos en formato wide
data_wide %>%
dplyr::slice(1:4) %>%
knitr::kable(., digits = 4)
id_j i1 i2 i3 i4 i5
1 -0.6496 -0.5137 -0.6978 0.0091 0.3443
2 1.8108 0.3455 0.8190 0.9825 0.9296
3 -1.3657 -1.0181 -1.3612 -1.0542 -0.5631
4 -0.5559 -0.3553 -0.6552 0.0533 -0.0303

3. Ajustar modelos

3.1 Modelo MLM tradicional

Empleamos la función lme4::lmer() para ajustar un modelo multinivel de tipo nulo, sobre los datos data_long.rds.

# ----------------------------------------------- 
# method 1: lem4::lmer()
# -----------------------------------------------

fit_null <- lme4::lmer(y ~ 1 + (1|id_j), data = data_long, REML = FALSE)
summary(fit_null)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: y ~ 1 + (1 | id_j)
   Data: data_long

     AIC      BIC   logLik deviance df.resid 
  3315.2   3332.6  -1654.6   3309.2     2497 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0237 -0.6097 -0.0063  0.6006  3.7461 

Random effects:
 Groups   Name        Variance Std.Dev.
 id_j     (Intercept) 0.9564   0.9780  
 Residual             0.1013   0.3183  
Number of obs: 2500, groups:  id_j, 500

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.006864   0.044196   0.155

3.2 Modelo MLM en lavaan

Empleamos la función lavaan::sem() para ajustar un modelo multinivel de tipo nulo, sobre los datos data_long.rds.

# ----------------------------------------------- 
# method 2: lavaan::sem()
# -----------------------------------------------

lavaan_mlm <- '
level: 1
# variance w
y~~w*y

level: 2
# variance b
y~~b*y

y~mu*1

# model terms
y00   := mu
var_w := w
var_b := b
icc   := b/(b+w)

'

fit_mlm <- lavaan::sem(model = lavaan_mlm, 
           data = data_long, 
           cluster = "id_j"
           )

summary(fit_mlm)
Length  Class   Mode 
     1 lavaan     S4 
lavaan::parameterEstimates(fit_mlm) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
dplyr::select(label, est, se) %>%
knitr::kable(., digits = 3)
label est se
y00 0.007 0.044
var_w 0.101 0.003
var_b 0.956 0.062
icc 0.904 0.006

3.2 Modelo MLM por medio de un modelo CFA

Empleamos la función lavaan::sem() para ajustar un modelo multinivel de tipo nulo a través de un modelo CFA, sobre los datos data_wide.

# ----------------------------------------------- 
# method 3: lavaan::sem()
# -----------------------------------------------

lavaan_cfa <- '
# random term
eta =~ 1*i1 
eta =~ 1*i2 
eta =~ 1*i3 
eta =~ 1*i4 
eta =~ 1*i5 

# latent mean
eta ~ mu*1

# person intercepts
i1  ~ d1 *1
i2  ~ d2 *1
i3  ~ d3 *1
i4  ~ d4 *1
i5  ~ d5 *1

# variance
eta~~b*eta

# residual variance
i1 ~~ w*i1
i2 ~~ w*i2
i3 ~~ w*i3
i4 ~~ w*i4
i5 ~~ w*i5

# model constraints

cw := (d1 +d2 +d3 +d4 +d5)/5

cw == 0

# model terms
y00   := mu
var_w := w
var_b := b
icc   := b/(b+w)


'

fit_cfa <- lavaan::sem(model = lavaan_cfa, 
           data = data_wide
           )

summary(fit_cfa)
Length  Class   Mode 
     1 lavaan     S4 
lavaan::parameterEstimates(fit_cfa) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
dplyr::select(label, est, se) %>%
knitr::kable(., digits = 3)
label est se
y00 0.007 0.044
var_w 0.101 0.003
var_b 0.956 0.062
icc 0.905 0.006

4. Tabla comparada de resultados

4.1 Extraer los estimados de cada modelo ajustado

# ----------------------------------------------- 
# retrive estimates
# -----------------------------------------------

# function to adjust to two decimal places
decimal <- function (x, k) {
    format(round(x, k), nsmall = k)
}

# lmer4::lmer()
table_est_1 <- broom.mixed::tidy(fit_null) %>%
mutate(label = c('y00','var_b', 'var_w')) %>%
mutate(est = case_when(
effect == 'fixed'    ~ estimate,
effect == 'ran_pars' ~ estimate*estimate
)) %>%
mutate(est_1_se = paste0(decimal(est, 2), ' (', decimal(std.error,2), ')')) %>%
dplyr::select(label, est_1_se)

# lavaan mlm
table_est_2 <- lavaan::parameterEstimates(fit_mlm) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
mutate(est_2_se = paste0(decimal(est, 2), ' (', decimal(se,2), ')')) %>%
dplyr::select(label, est_2_se) 

# lavaan cfa
table_est_3 <- lavaan::parameterEstimates(fit_cfa) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
mutate(est_3_se = paste0(decimal(est, 2), ' (', decimal(se,2), ')')) %>%
dplyr::select(label, est_3_se) 

4.2 Tabla comparada

# ----------------------------------------------- 
# display table comparison
# -----------------------------------------------

dplyr::full_join(table_est_1, table_est_2, by = 'label') %>%
dplyr::full_join(., table_est_3, by = 'label') %>%
rename(
`lme4::lmer()` = 2,
`lavaan::sem()-MLM` = 3,
`lavaan::sem()-CFA` = 4
) %>%
knitr::kable()
label lme4::lmer() lavaan::sem()-MLM lavaan::sem()-CFA
y00 0.01 (0.04) 0.01 (0.04) 0.01 (0.04)
var_b 0.96 ( NA) 0.96 (0.06) 0.96 (0.06)
var_w 0.10 ( NA) 0.10 (0.00) 0.10 (0.00)
icc 0.90 (0.01) 0.90 (0.01)

Anexos

Código Mplus para generar los datos simulados

TITLE:
null mlm;

MONTECARLO:
 
NAMES         = y;       
NOBSERVATIONS = 2500;    
NCSIZES       = 1;       
CSIZES        = 500 (5); 
NREPS         = 1;       
SEED          = 4046;    
SAVE          = icc_90.dat;

MODEL POPULATION:
 
%WITHIN%                 
!residual variance       
y*0.1;

%BETWEEN%               
!latent mean            
[y@0];                  
y*0.9;

ANALYSIS:

TYPE = TWOLEVEL;
ESTIMATOR = ML;

MODEL:
 
%WITHIN%                 
!residual variance       
y*0.1;

%BETWEEN%               
!latent mean            
[y@0];                  
y*0.9;

OUTPUT:
TECH9;